# Analyzing metabolomics data for Mutographs ----
# 02_regressions
# Running regressions models

rm(list = ls())

# 0 - Definition of libraries, paths and functions ----

Sys.setlocale("LC_TIME", "C")

library(stringr)
library(tidyr)
library(Hmisc)
library(MASS)
library(data.table)
library(dplyr)
library(forcats)
library(tibble)
library(openxlsx)
library(ggplot2)
library(cowplot)
library(scales)
library(patchwork)
library(ggpubr)

# Loading previously processed data
metabo_RAW <- read.csv("data/metabolomics_normalized_data.csv", check.names=FALSE)
metabo_PROC <- readRDS("output/metabo_PROC")
metabo_MAN <- readRDS("output/metabo_MAN")
metabo_pearson <- readRDS("output/metabo_pearson")
proxy_table <- readRDS("output/proxy_table")

features <- names(metabo_RAW)[str_detect(names(metabo_RAW),"@")]
metabolites <- names(metabo_PROC)[str_detect(names(metabo_PROC),"@")]

cosmic_attrib <- read.csv("data/sigs/pruned_attribution_RCC_Manuscript_COSMIC_SBS96_abs_mutations.csv") %>% 
  rename(donor_id=X, SBS1536A_cosmic = SBS1536A, SBS1536B_cosmic = SBS1536B, SBS1536F_cosmic = SBS1536F, SBS1536I_cosmic = SBS1536I)
dbs_attrib <- read.csv("data/sigs/pruned_attribution_RCC_Manuscript_COSMIC_DBS78_abs_mutations.csv") %>% rename(donor_id=X)
id_attrib <- read.csv("data/sigs/pruned_attribution_RCC_Manuscript_COSMIC_ID83_abs_mutations.csv") %>% rename(donor_id=X)
cn_attrib <- read.csv("data/sigs/pruned_attribution_RCC_Manuscript_COSMIC_CNV48_abs_mutations.csv") %>% rename(donor_id=X)

sbs_dn_attrib <- read.csv("data/sigs/pruned_attribution_RCC_Manuscript_denovo_SBS1536_abs_mutations.csv") %>% rename(donor_id=X)
dbs_dn_attrib <- read.csv("data/sigs/pruned_attribution_RCC_Manuscript_denovo_DBS78_abs_mutations.csv") %>% rename(donor_id=X)
id_dn_attrib <- read.csv("data/sigs/pruned_attribution_RCC_Manuscript_denovo_ID83_abs_mutations.csv") %>% rename(donor_id=X)
cn_dn_attrib <- read.csv("data/sigs/pruned_attribution_RCC_Manuscript_denovo_CNV48_abs_mutations.csv") %>% rename(donor_id=X)
sv_dn_attrib <- read.csv("data/sigs/pruned_attribution_RCC_Manuscript_denovo_SV32_abs_mutations.csv") %>% rename(donor_id=X)

sigs_cosmic <- names(cosmic_attrib)[2:15]
dbs_cosmic <- names(dbs_attrib)[2:6]
id_cosmic <- names(id_attrib)[2:10]
cn_cosmic <- names(cn_attrib)[2:5]

sigs_dn <- names(sbs_dn_attrib)[2:13]
dbs_dn <- names(dbs_dn_attrib)[2:5]
id_dn <- names(id_dn_attrib)[2:8]
cn_dn <- names(cn_dn_attrib)[2:4]
sv_dn <- names(sv_dn_attrib)[2:4]

# Loading data from curated TMAP peak
tmap_corrected <- read.csv("data/tmap_corrected_peak.csv")

# Checking AA and SBS12 exposed cases
# Listing cases with >10% of signature attribution due to SBS22a + SBS22b or SBS12
cosmic_relative <- cosmic_attrib %>% 
  mutate(total_attrib = rowSums(.[2:15])) %>% 
  mutate(across(starts_with("SBS"), ~./total_attrib))

aa_cases <- filter(cosmic_relative, SBS22 + SBS1536I_cosmic > 0.1)$donor_id
SBS12_cases <- filter(cosmic_relative, SBS12 > 0.1)$donor_id  

# Running regressions - All cases ----
# * COSMIC ----
# ** SBS ----
# With quasipoisson or logistic regression / COSMIC signatures
sparseness_sigcosm <- colSums(select(metabo_MAN, all_of(sigs_cosmic))==0)/NROW(metabo_MAN)*100
reg_cosm <- list()
for (sig in sigs_cosmic){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    
    # Model: Quasipoisson if sparseness below 30%, logistic (A/B median) else
    if (sparseness_sigcosm[sig] >= 30){
      lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    } else {
      lreg <- glm(get(sig) ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "quasipoisson")
    }
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_cosm[[sig]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=sig, pval_adj = p.adjust(pval, "bonferroni"))
}
# Adding logistic for SBS1 and SBS1536A/B (cosmic) as sparseness is < 30%. Others are too zero-inflated (>50%)
for (sig in c("SBS1", "SBS1536A_cosmic", "SBS1536B_cosmic")){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_cosm[[str_c(sig, "_logistic")]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=str_c(sig, "_logistic"), pval_adj = p.adjust(pval, "bonferroni")) 
}
# saveRDS(reg_cosm, "results/regressions/reg_cosm")

res_reg_cosm <- rbindlist(reg_cosm) %>% 
  group_by(signature) %>% 
  mutate(sig_label = factor(case_when(signature == "SBS1536B_cosmic" ~ "SBS40a", signature == "SBS1536A_cosmic" ~ "SBS40b", signature == "SBS1536F_cosmic" ~ "SBS40c",
                                      signature == "SBS22" ~ "SBS22a", signature == "SBS1536I_cosmic" ~ "SBS22b", TRUE ~ signature),
                            levels = c("SBS1","SBS2","SBS4","SBS5","SBS12","SBS13","SBS18","SBS21","SBS22a","SBS22b","SBS40a","SBS40b","SBS40c","SBS44")),
         label = ifelse(pval < 0.05/length(metabolites), str_c(metabo, sig_label, sep = " \n "), ""),
         diffexpressed = case_when(pval < 0.05/length(metabolites) & Estimate > 0 ~ 1,
                                   pval < 0.05/length(metabolites) & Estimate <= 0 ~ -1,
                                   TRUE ~ 0))
# saveRDS(res_reg_cosm, "results/regressions/res_reg_cosm")

# ** DBS ----
sparseness_dbscosm <- colSums(select(filter(metabo_MAN, !is.na(DBS2)), all_of(dbs_cosmic))==0)/NROW(filter(metabo_MAN, !is.na(DBS2)))*100
# All are too sparse for quasipoisson (> 30% 0s)
reg_dbs_cosm <- list()
for (sig in dbs_cosmic){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN %>% filter(!is.na(DBS2)) %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    
    # Model: Quasipoisson if sparseness below 30%, logistic (A/B median) else
    if (sparseness_dbscosm[sig] >= 30){
      lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    } else {
      lreg <- glm(get(sig) ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "quasipoisson")
    }
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_dbs_cosm[[sig]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=sig, pval_adj = p.adjust(pval, "bonferroni"))
}
# saveRDS(reg_dbs_cosm, "results/regressions/reg_dbs_cosm")
res_reg_dbs_cosm <- rbindlist(reg_dbs_cosm) %>% 
  group_by(signature) %>% 
  mutate(label = ifelse(pval < 0.05/length(metabolites), str_c(metabo, signature, sep = " \n "), ""),
         diffexpressed = case_when(pval < 0.05/length(metabolites) & Estimate > 0 ~ 1,
                                   pval < 0.05/length(metabolites) & Estimate <= 0 ~ -1,
                                   TRUE ~ 0))
# saveRDS(res_reg_dbs_cosm, "results/regressions/res_reg_dbs_cosm")

# ** ID ----
sparseness_idcosm <- colSums(select(metabo_MAN, all_of(id_cosmic))==0)/NROW(metabo_MAN)*100
# All are sparse except ID1 (14% 0s), ID5 (6% 0s) and ID8 (27% 0s)
reg_id_cosm <- list()
for (sig in id_cosmic){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    
    # Model: Quasipoisson if sparseness below 30%, logistic (A/B median) else
    if (sparseness_idcosm[sig] >= 30){
      lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    } else {
      lreg <- glm(get(sig) ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "quasipoisson")
    }
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_id_cosm[[sig]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=sig, pval_adj = p.adjust(pval, "bonferroni"))
}
# Adding logistic for ID1, ID5 and ID8 as sparseness is < 30%. Others are too zero-inflated (>50%)
for (sig in c("ID1", "ID5", "ID8")){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_id_cosm[[str_c(sig, "_logistic")]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=str_c(sig, "_logistic"), pval_adj = p.adjust(pval, "bonferroni")) 
}
# saveRDS(reg_id_cosm, "results/regressions/reg_id_cosm")
res_reg_id_cosm <- rbindlist(reg_id_cosm) %>% 
  group_by(signature) %>% 
  mutate(label = ifelse(pval < 0.05/length(metabolites), str_c(metabo, signature, sep = " \n "), ""),
         diffexpressed = case_when(pval < 0.05/length(metabolites) & Estimate > 0 ~ 1,
                                   pval < 0.05/length(metabolites) & Estimate <= 0 ~ -1,
                                   TRUE ~ 0))
# saveRDS(res_reg_id_cosm, "results/regressions/res_reg_id_cosm")

# ** CNV ----
metabo_MAN_CN <- metabo_MAN %>% left_join(cn_attrib) %>% left_join(cn_dn_attrib) %>% # Adding CNV COSMIC/DN signatures
  filter(!is.na(CN1)) %>% # Removing all cases for which attribution was not performed
  mutate(across(
    any_of(c(cn_cosmic, cn_dn, sv_dn)),
    .fns = list(cat = ~cut2(.,g=2), int = ~qnorm((rank(.,na.last="keep")-0.5)/sum(!is.na(.))), logdelta = ~log2(.+1)),
    .names = "{.col}_{.fn}"))
sparseness_cncosm <- colSums(select(metabo_MAN_CN, all_of(cn_cosmic))==0)/NROW(metabo_MAN_CN)*100
# All are sparse except CN1 (15% 0s)
reg_cn_cosm <- list()
for (sig in cn_cosmic){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN_CN %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    
    # Model: Quasipoisson if sparseness below 30%, logistic (A/B median) else
    if (sparseness_cncosm[sig] >= 30){
      lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    } else {
      lreg <- glm(get(sig) ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "quasipoisson")
    }
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_cn_cosm[[sig]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=sig, pval_adj = p.adjust(pval, "bonferroni"))
}
# Adding logistic for CN1 as sparseness is < 30%. Others are too zero-inflated (>50%)
for (sig in c("CN1")){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN_CN %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_cn_cosm[[str_c(sig, "_logistic")]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=str_c(sig, "_logistic"), pval_adj = p.adjust(pval, "bonferroni")) 
}
# saveRDS(reg_cn_cosm, "results/regressions/reg_cn_cosm")
res_reg_cn_cosm <- rbindlist(reg_cn_cosm) %>% 
  group_by(signature) %>% 
  mutate(label = ifelse(pval < 0.05/length(metabolites), str_c(metabo, signature, sep = " \n "), ""),
         diffexpressed = case_when(pval < 0.05/length(metabolites) & Estimate > 0 ~ 1,
                                   pval < 0.05/length(metabolites) & Estimate <= 0 ~ -1,
                                   TRUE ~ 0))
# saveRDS(res_reg_cn_cosm, "results/regressions/res_reg_cn_cosm")

# ** Mutational burden ----
mutburden_sbs <- read.csv("data/sigs/output_RCC_Manuscript_COSMIC_SBS96_stat_info.csv") %>% 
  rename(donor_id=Sample, mutburden_sbs=Mutational.burden) %>% select(donor_id, mutburden_sbs)
mutburden_dbs <- read.csv("data/sigs/output_RCC_Manuscript_COSMIC_DBS78_stat_info.csv") %>% 
  rename(donor_id=Sample, mutburden_dbs=Mutational.burden) %>% select(donor_id, mutburden_dbs)
mutburden_id <- read.csv("data/sigs/output_RCC_Manuscript_COSMIC_ID83_stat_info.csv") %>% 
  rename(donor_id=Sample, mutburden_id=Mutational.burden) %>% select(donor_id, mutburden_id)
mutburden_cn <- read.csv("data/sigs/output_RCC_Manuscript_COSMIC_CNV48_stat_info.csv") %>% 
  rename(donor_id=Sample, mutburden_cn=Mutational.burden) %>% select(donor_id, mutburden_cn)
mutburden_sv <- read.csv("data/sigs/output_RCC_Manuscript_denovo_SV32_stat_info.csv") %>% 
  rename(donor_id=Sample, mutburden_sv=Mutational.burden) %>% select(donor_id, mutburden_sv)
mutburden <- mutburden_sbs %>% left_join(mutburden_dbs) %>% left_join(mutburden_id) %>% left_join(mutburden_cn) %>% left_join(mutburden_sv)

mutburden_all <- names(mutburden)[2:6]

metabo_MAN_burden <- metabo_MAN %>% left_join(mutburden)
# Using only Quasipoisson regression for burden as missing for none
reg_burden <- list()
for (sig in mutburden_all){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN_burden %>% filter(!is.na(get(sig))) %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    
    lreg <- glm(get(sig) ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "quasipoisson")
    
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_burden[[sig]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=sig, pval_adj = p.adjust(pval, "bonferroni"))
}
# saveRDS(reg_burden, "results/regressions/reg_burden")
res_reg_burden<- rbindlist(reg_burden) %>% 
  group_by(signature) %>% 
  mutate(label = ifelse(pval < 0.05/length(metabolites), str_c(metabo, signature, sep = " \n "), ""),
         diffexpressed = case_when(pval < 0.05/length(metabolites) & Estimate > 0 ~ 1,
                                   pval < 0.05/length(metabolites) & Estimate <= 0 ~ -1,
                                   TRUE ~ 0))
# saveRDS(res_reg_burden, "results/regressions/res_reg_burden")

# ** Mutational burden - Excluding Romania and Serbia ----
# Using only Quasipoisson regression for burden as missing for none
reg_burden_rsexcluded <- list()
for (sig in mutburden_all){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN_burden %>% filter(!is.na(get(sig)) & country %nin% c("Romania","Serbia")) %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    
    lreg <- glm(get(sig) ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "quasipoisson")
    
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_burden_rsexcluded[[sig]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=sig, pval_adj = p.adjust(pval, "bonferroni"))
}
# saveRDS(reg_burden_rsexcluded, "results/regressions/reg_burden_rsexcluded")
res_reg_burden_rsexcluded<- rbindlist(reg_burden_rsexcluded) %>% 
  group_by(signature) %>% 
  mutate(label = ifelse(pval < 0.05/length(metabolites), str_c(metabo, signature, sep = " \n "), ""),
         diffexpressed = case_when(pval < 0.05/length(metabolites) & Estimate > 0 ~ 1,
                                   pval < 0.05/length(metabolites) & Estimate <= 0 ~ -1,
                                   TRUE ~ 0))
# saveRDS(res_reg_burden_rsexcluded, "results/regressions/res_reg_burden_rsexcluded")

# ** Mutational burden - Excluding AA and SBS12 > 10% cases ----
# Using only Quasipoisson regression for burden as missing for none
reg_burden_aa12excluded <- list()
for (sig in mutburden_all){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN_burden %>% filter(!is.na(get(sig)) & !donor_id %in% c(aa_cases,SBS12_cases)) %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    
    lreg <- glm(get(sig) ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "quasipoisson")
    
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_burden_aa12excluded[[sig]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=sig, pval_adj = p.adjust(pval, "bonferroni"))
}
# saveRDS(reg_burden_aa12excluded, "results/regressions/reg_burden_aa12excluded")
res_reg_burden_aa12excluded<- rbindlist(reg_burden_aa12excluded) %>% 
  group_by(signature) %>% 
  mutate(label = ifelse(pval < 0.05/length(metabolites), str_c(metabo, signature, sep = " \n "), ""),
         diffexpressed = case_when(pval < 0.05/length(metabolites) & Estimate > 0 ~ 1,
                                   pval < 0.05/length(metabolites) & Estimate <= 0 ~ -1,
                                   TRUE ~ 0))
# saveRDS(res_reg_burden_aa12excluded, "results/regressions/res_reg_burden_aa12excluded")

# * De novo ----
# ** SBS ----
sparseness_sigdn <- colSums(select(metabo_MAN, all_of(sigs_dn))==0)/NROW(metabo_MAN)*100
# All signatures are too sparse for quasipoisson (> 30% 0s)
reg_dn <- list()
for (sig in sigs_dn){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    
    # Model: Quasipoisson if sparseness below 30%, logistic (A/B median) else
    if (sparseness_sigdn[sig] >= 30){
      lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    } else {
      lreg <- glm(get(sig) ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "quasipoisson")
    }
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_dn[[sig]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=sig, pval_adj = p.adjust(pval, "bonferroni"))
}
# saveRDS(reg_dn, "results/regressions/reg_dn")

res_reg_dn<- rbindlist(reg_dn) %>% 
  group_by(signature) %>% 
  mutate(label = ifelse(pval < 0.05/length(metabolites), str_c(metabo, signature, sep = " \n "), ""),
         diffexpressed = case_when(pval < 0.05/length(metabolites) & Estimate > 0 ~ 1,
                                   pval < 0.05/length(metabolites) & Estimate <= 0 ~ -1,
                                   TRUE ~ 0))
# saveRDS(res_reg_dn, "results/regressions/res_reg_dn")

# ** DBS ----
# Adding de novo DBS and ID in the compiled dataset
metabo_MAN_DN <- metabo_MAN %>% select(-DBS78C, -DBS78D, -ID83C) %>% left_join(dbs_dn_attrib) %>% left_join(id_dn_attrib) %>% # Removing COSMIC versions of de novo signatures
  mutate(across(
    any_of(c(dbs_dn,id_dn)),
    .fns = list(cat = ~cut2(.,g=2), int = ~qnorm((rank(.,na.last="keep")-0.5)/sum(!is.na(.))), logdelta = ~log2(.+1)),
    .names = "{.col}_{.fn}"))
# With DBS de novo sigs
sparseness_dbsdn <- colSums(select(filter(metabo_MAN_DN, !is.na(DBS78A)), all_of(dbs_dn))==0)/NROW(filter(metabo_MAN_DN, !is.na(DBS78A)))*100
# All are too sparse for quasipoisson (> 30% 0s)
reg_dbs_dn <- list()
for (sig in dbs_dn){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN_DN %>% filter(!is.na(DBS78A)) %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    
    # Model: Quasipoisson if sparseness below 30%, logistic (A/B median) else
    if (sparseness_dbsdn[sig] >= 30){
      lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    } else {
      lreg <- glm(get(sig) ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "quasipoisson")
    }
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_dbs_dn[[sig]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=sig, pval_adj = p.adjust(pval, "bonferroni"))
}
# saveRDS(reg_dbs_dn, "results/regressions/reg_dbs_dn")
res_reg_dbs_dn <- rbindlist(reg_dbs_dn) %>% 
  group_by(signature) %>% 
  mutate(label = ifelse(pval < 0.05/length(metabolites), str_c(metabo, signature, sep = " \n "), ""),
         diffexpressed = case_when(pval < 0.05/length(metabolites) & Estimate > 0 ~ 1,
                                   pval < 0.05/length(metabolites) & Estimate <= 0 ~ -1,
                                   TRUE ~ 0))
# saveRDS(res_reg_dbs_dn, "results/regressions/res_reg_dbs_dn")

# ** ID ----
sparseness_iddn <- colSums(select(metabo_MAN_DN, all_of(id_dn))==0)/NROW(metabo_MAN_DN)*100
# All are sparse except ID83A (24% 0s) and ID83F (23% 0s)
reg_id_dn <- list()
for (sig in id_dn){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN_DN %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    
    # Model: Quasipoisson if sparseness below 30%, logistic (A/B median) else
    if (sparseness_iddn[sig] >= 30){
      lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    } else {
      lreg <- glm(get(sig) ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "quasipoisson")
    }
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_id_dn[[sig]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=sig, pval_adj = p.adjust(pval, "bonferroni"))
}
# Adding logistic for ID83A and ID83F as sparseness is < 30% (to have both quasipoisson and logistic)
for (sig in c("ID83A", "ID83F")){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN_DN %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_id_dn[[str_c(sig, "_logistic")]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=str_c(sig, "_logistic"), pval_adj = p.adjust(pval, "bonferroni")) 
}
# saveRDS(reg_id_dn, "results/regressions/reg_id_dn")
res_reg_id_dn <- rbindlist(reg_id_dn) %>%
  group_by(signature) %>% 
  mutate(label = ifelse(pval < 0.05/length(metabolites), str_c(metabo, signature, sep = " \n "), ""),
         diffexpressed = case_when(pval < 0.05/length(metabolites) & Estimate > 0 ~ 1,
                                   pval < 0.05/length(metabolites) & Estimate <= 0 ~ -1,
                                   TRUE ~ 0))
# saveRDS(res_reg_id_dn, "results/regressions/res_reg_id_dn")

# ** CNV ----
sparseness_cndn <- colSums(select(metabo_MAN_CN, all_of(cn_dn))==0)/NROW(metabo_MAN_CN)*100
# All are too sparse for quasipoisson except SBSCNVA (12.25% 0s). Rest is > 30% 0s.
reg_cn_dn <- list()
for (sig in cn_dn){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN_CN %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    
    # Model: Quasipoisson if sparseness below 30%, logistic (A/B median) else
    if (sparseness_cndn[sig] >= 30){
      lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    } else {
      lreg <- glm(get(sig) ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "quasipoisson")
    }
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_cn_dn[[sig]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=sig, pval_adj = p.adjust(pval, "bonferroni"))
}
# Adding logistic for SBSCNVA as sparseness is < 30%. Others are too zero-inflated (>50%)
for (sig in c("SBSCNVA")){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN_CN %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_cn_dn[[str_c(sig, "_logistic")]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=str_c(sig, "_logistic"), pval_adj = p.adjust(pval, "bonferroni")) 
}
# saveRDS(reg_cn_dn, "results/regressions/reg_cn_dn")
res_reg_cn_dn <- rbindlist(reg_cn_dn) %>% 
  group_by(signature) %>% 
  mutate(label = ifelse(pval < 0.05/length(metabolites), str_c(metabo, signature, sep = " \n "), ""),
         diffexpressed = case_when(pval < 0.05/length(metabolites) & Estimate > 0 ~ 1,
                                   pval < 0.05/length(metabolites) & Estimate <= 0 ~ -1,
                                   TRUE ~ 0))
# saveRDS(res_reg_cn_dn, "results/regressions/res_reg_cn_dn")

# ** SV ----
metabo_MAN_SV <- metabo_MAN %>% left_join(sv_dn_attrib) %>% # Adding SV DN signatures
  mutate(across(
    any_of(c(sv_dn)),
    .fns = list(cat = ~cut2(.,g=2), int = ~qnorm((rank(.,na.last="keep")-0.5)/sum(!is.na(.))), logdelta = ~log2(.+1)),
    .names = "{.col}_{.fn}"))
sparseness_svdn <- colSums(select(metabo_MAN_SV, all_of(sv_dn))==0)/NROW(metabo_MAN_SV)*100
# All are too sparse for quasipoisson except SBSSVA (14.2% 0s). Rest is > 30% 0s.
reg_sv_dn <- list()
for (sig in sv_dn){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN_SV %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    
    # Model: Quasipoisson if sparseness below 30%, logistic (A/B median) else
    if (sparseness_svdn[sig] >= 30){
      lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    } else {
      lreg <- glm(get(sig) ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "quasipoisson")
    }
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_sv_dn[[sig]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=sig, pval_adj = p.adjust(pval, "bonferroni"))
}
# Adding logistic for SBSSVA as sparseness is < 30%. Others are too zero-inflated (>50%)
for (sig in c("SBSSVA")){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN_SV %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_sv_dn[[str_c(sig, "_logistic")]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=str_c(sig, "_logistic"), pval_adj = p.adjust(pval, "bonferroni")) 
}
# saveRDS(reg_sv_dn, "results/regressions/reg_sv_dn")
res_reg_sv_dn <- rbindlist(reg_sv_dn) %>% 
  group_by(signature) %>% 
  mutate(label = ifelse(pval < 0.05/length(metabolites), str_c(metabo, signature, sep = " \n "), ""),
         diffexpressed = case_when(pval < 0.05/length(metabolites) & Estimate > 0 ~ 1,
                                   pval < 0.05/length(metabolites) & Estimate <= 0 ~ -1,
                                   TRUE ~ 0))
# saveRDS(res_reg_sv_dn, "results/regressions/res_reg_sv_dn")

# Running regressions - Only in Serbia & Romania for SBS22 and SBS1536I (SBS22b) ----

sparseness_sigcosm_ro <- colSums(select(filter(metabo_MAN, country %in% c("Romania", "Serbia")), all_of(sigs_cosmic))==0)/NROW(filter(metabo_MAN, country %in% c("Romania", "Serbia")))*100
# SBS22 and SBS1536I still super sparse in Romania and Serbia (> 40%), so logistic regression is used
reg_cosm_sro <- list()
for (sig in c("SBS22", "SBS1536I_cosmic")){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN %>% filter(country %in% c("Romania", "Serbia")) %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    
    # Model: Quasipoisson if sparseness below 30%, logistic (P/A) else
    if (sparseness_sigcosm_ro[sig] >= 30){
      # print("Logistic")
      lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    } else {
      # print("Quasipoisson")
      lreg <- glm(get(sig) ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "quasipoisson")
    }
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_cosm_sro[[sig]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=sig, pval_adj = p.adjust(pval, "bonferroni"))
}
# saveRDS(reg_cosm_sro, "results/regressions/reg_cosm_sro")

res_reg_cosm_sro <- rbindlist(reg_cosm_sro) %>% 
  group_by(signature) %>% 
  mutate(sig_label = factor(case_when(signature == "SBS1536B_cosmic" ~ "SBS40a", signature == "SBS1536A_cosmic" ~ "SBS40b", signature == "SBS1536F_cosmic" ~ "SBS40c",
                                      signature == "SBS22" ~ "SBS22a", signature == "SBS1536I_cosmic" ~ "SBS22b", TRUE ~ signature),
                            levels = c("SBS1","SBS2","SBS4","SBS5","SBS12","SBS13","SBS18","SBS21","SBS22a","SBS22b","SBS40a","SBS40b","SBS40c","SBS44")),
         label = ifelse(pval < 0.05/length(metabolites), str_c(metabo, sig_label, sep = " \n "), ""),
         diffexpressed = case_when(pval < 0.05/length(metabolites) & Estimate > 0 ~ 1,
                                   pval < 0.05/length(metabolites) & Estimate <= 0 ~ -1,
                                   TRUE ~ 0))
# saveRDS(res_reg_cosm_sro, "results/regressions/res_reg_cosm_sro")

# Same for DBS78D and ID83C
sparseness_dbscosm_ro <- colSums(select(filter(metabo_MAN, country %in% c("Romania", "Serbia") & !is.na(DBS2)), DBS78D)==0)/NROW(filter(metabo_MAN, country %in% c("Romania", "Serbia") & !is.na(DBS2)))*100
sparseness_idcosm_ro <- colSums(select(filter(metabo_MAN, country %in% c("Romania", "Serbia")), ID83C)==0)/NROW(filter(metabo_MAN, country %in% c("Romania", "Serbia")))*100

# DBS78D: 67% of 0s / ID83C: 89% of 0s 
# Both are too sparse for quasipoisson (> 30% 0s) so using logistic

# DBS78D
reg_dbs_cosm_ro <- list()
for (sig in c("DBS78D")){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN %>% filter(country %in% c("Romania", "Serbia") & !is.na(DBS2)) %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    
    # Model: Logistic (A/B median)
    lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_dbs_cosm_ro[[sig]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=sig, pval_adj = p.adjust(pval, "bonferroni"))
}
# saveRDS(reg_dbs_cosm_ro, "results/regressions/reg_dbs_cosm_ro")

res_reg_dbs_cosm_sro <- rbindlist(reg_dbs_cosm_ro) %>% 
  group_by(signature) %>% 
  mutate(label = ifelse(pval < 0.05/length(metabolites), str_c(metabo, sig_label, sep = " \n "), ""),
         diffexpressed = case_when(pval < 0.05/length(metabolites) & Estimate > 0 ~ 1,
                                   pval < 0.05/length(metabolites) & Estimate <= 0 ~ -1,
                                   TRUE ~ 0))
# saveRDS(res_reg_dbs_cosm_sro, "results/regressions/res_reg_dbs_cosm_sro")
# Nothing significant after adjustment

# ID83C
reg_id_cosm_ro <- list()
for (sig in c("ID83C")){
  ltmp <- list()
  for (metabolite in metabolites){
    
    df_mod <- metabo_MAN %>% filter(country %in% c("Romania", "Serbia")) %>% select(metabolite, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    
    # Model: Logistic (A/B median)
    lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = metabolite)
    ltmp[[metabolite]] <- resreg
  }
  reg_id_cosm_ro[[sig]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=sig, pval_adj = p.adjust(pval, "bonferroni"))
}
# saveRDS(reg_id_cosm_ro, "results/regressions/reg_id_cosm_ro")

res_reg_id_cosm_sro <- rbindlist(reg_id_cosm_ro) %>% 
  group_by(signature) %>% 
  mutate(label = ifelse(pval < 0.05/length(metabolites), str_c(metabo, sig_label, sep = " \n "), ""),
         diffexpressed = case_when(pval < 0.05/length(metabolites) & Estimate > 0 ~ 1,
                                   pval < 0.05/length(metabolites) & Estimate <= 0 ~ -1,
                                   TRUE ~ 0))
# saveRDS(res_reg_id_cosm_sro, "results/regressions/res_reg_id_cosm_sro")
# Nothing significant after adjustment


# Running regression with the curated TMAP peak, as well as with Creatitine ----

# TMAP corrected manually by Pekka
tmap_processed <- tmap_corrected %>% 
  # Negative to NA
  mutate(across(where(is.numeric), ~ifelse(.<0,NA,.))) %>%
  # Missing to 1/5 of minimum
  mutate(across(where(is.numeric), ~replace_na(.,1/5*min(.,na.rm=T)))) %>%
  # Log transformation
  mutate(across(where(is.numeric), ~log10(.))) %>% 
  # Pareto scaling
  mutate(across(where(is.numeric), ~(. - mean(.,na.rm=T))/sqrt(sd(.,na.rm=T))))

TMAP_RAW <- metabo_RAW %>% 
  rename(peak240 = `240.1468@0.8929933`, creatinine=`113.0589@0.6194511`,
         cotinine = `176.095@1.5373601`, hydroxycotinine = `192.0896@0.8945959`,
         TMAP1 = `228.1464@0.85400295`, TMAP2 = `228.1466@0.85400516`, TMAP3 = `228.1463@0.8544947`) %>% 
  mutate(batch = relevel(factor(batch), ref="1")) %>% 
  left_join(tmap_corrected)

TMAP_MAN <- metabo_MAN %>% 
  rename(peak240 = `240.1468@0.8929933`, creatinine=`113.0589@0.6194511`,
         cotinine = `176.095@1.5373601`, hydroxycotinine = `192.0896@0.8945959`) %>% 
  left_join(tmap_processed) %>% 
  left_join(cn_attrib) %>% 
  left_join(mutburden)

reg_TMAP_240 <- glm(TMAP ~ peak240 + batch/acq_order, data = TMAP_RAW)
summary(reg_TMAP_240)$coefficients
reg_TMAP_CREA <- glm(TMAP ~ creatinine + batch/acq_order, data = TMAP_RAW)
summary(reg_TMAP_CREA)$coefficients

reg_TMAP <- glm(SBS1536A_cosmic ~ TMAP + age_diag + sex + bmi +  batch/acq_order, data = TMAP_MAN, family = "quasipoisson")
summary(reg_TMAP)
reg_CREA <- glm(SBS1536A_cosmic ~ creatinine + age_diag + sex + bmi +  batch/acq_order, data = TMAP_MAN, family = "quasipoisson")
summary(reg_CREA)

# Running regressions of TMAP (curated) vs all signatures (SBS/DBS/ID) ----
# Quasipoisson for SBS1, SBS1536A_cosmic, SBS1536B_cosmic, ID1, ID5 & ID8, CN1 - Rest is logistic
# Only in Romania & Serbia for SBS22, SBS1536I_cosmic, ID83C, DBS78D
# Removing AA & SBS12 cases for burdens
all_sigs <- c(sigs_cosmic, dbs_cosmic, id_cosmic, cn_cosmic, mutburden_all)
all_reg_TMAP <- list()
for(sig in all_sigs){
  # Selecting cases
  if (sig %in% c("SBS22", "SBS1536I_cosmic", "DBS78D", "ID83C")){
    df_mod <- filter(TMAP_MAN, country %in% c("Romania", "Serbia") & !is.na(get(sig)))
  } else {df_mod <- filter(TMAP_MAN, !is.na(get(sig)))}
  df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
  # Selecting model
  if (sig %in% c("SBS1", "SBS1536A_cosmic","SBS1536B_cosmic","ID1","ID5","ID8","CN1",mutburden_all)){
    lreg <- glm(get(sig) ~ TMAP + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "quasipoisson")
  } else {lreg <- glm(sig_cat ~ TMAP + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")}
  
  resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
  names(resreg) <- c("rowname","Estimate","pval")
  all_reg_TMAP[[sig]]<- resreg %>% filter(rowname == "TMAP") %>% mutate(pval_adj = min(1,pval*length(all_sigs)), signature = sig,
                                                                        model = case_when(signature %in% c("SBS22", "SBS1536I_cosmic", "DBS78D", "ID83C") ~ "Logistic - Romania & Serbia only",
                                                                                          signature %in% c("SBS1", "SBS1536A_cosmic","SBS1536B_cosmic","ID1","ID5","ID8","CN1",mutburden_all) ~ "Quasipoisson",
                                                                                          TRUE ~ "Logistic"))
}
res_all_reg_TMAP <- rbindlist(all_reg_TMAP)
# saveRDS(res_all_reg_TMAP, "results/regressions/res_all_reg_TMAP")

# Rerunning regressions for TMAP vs. burdens removing AA and SBS12 cases
all_reg_TMAP_AA12ex <- list()
for(sig in mutburden_all){
  # Selecting cases
  df_mod <- filter(TMAP_MAN, !is.na(get(sig)) & !donor_id %in% c(aa_cases,SBS12_cases))
  df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
  # Selecting model
  lreg <- glm(get(sig) ~ TMAP + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "quasipoisson")
  resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
  names(resreg) <- c("rowname","Estimate","pval")
  all_reg_TMAP_AA12ex[[sig]]<- resreg %>% filter(rowname == "TMAP") %>% mutate(pval_adj = min(1,pval*length(mutburden_all)), signature = sig,
                                                                               model = "Quasipoisson - Without AA exposed")
}
res_all_reg_TMAP_AA12ex <- rbindlist(all_reg_TMAP_AA12ex)
# saveRDS(res_all_reg_TMAP_AA12ex, "results/regressions/res_all_reg_TMAP_AA12ex")

# Rerunning regressions for TMAP vs. burdens with all cases
TMAP_MAN_burden <- TMAP_MAN %>% 
  mutate(mutburden_sbs_m40b = mutburden_sbs - SBS1536A)
all_reg_TMAP_burden <- list()
for(sig in c(mutburden_all, "mutburden_sbs_m40b")){
  # Selecting cases
  df_mod <- TMAP_MAN_burden %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
  # Selecting model
  lreg <- glm(get(sig) ~ TMAP + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "quasipoisson")
  resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
  names(resreg) <- c("rowname","Estimate","pval")
  all_reg_TMAP_burden[[sig]]<- resreg %>% filter(rowname == "TMAP") %>% mutate(pval_adj = min(1,pval*length(mutburden_all)), signature = sig,
                                                                               model = "Quasipoisson")
}
res_all_reg_TMAP_burden <- rbindlist(all_reg_TMAP_burden)
# saveRDS(res_all_reg_TMAP_AA12ex, "results/regressions/all_reg_TMAP_burden")